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ABSTRACT 

We perform a stability analysis of a tidally excited nonlinear internal gravity wave 
near the centre of a solar-type star in two-dimensional cylindrical geometry. The mo- 
tivation is to understand the tidal interaction between short-period planets and their 
slowly rotating solar-type host stars, which involves the launching of internal gravity 
waves at the top of the radiation zone that propagate towards the centre of the star. 
Studying the instabilities of these waves near the centre, where nonlinearities are most 
important, is essential, since it may have implications for the survival of short-period 
planets orbiting solar- type stars. When these waves have sufficient amplitude to over- 
turn the stratification, they break and form a critical layer, which efficiently absorbs 
subsequent ingoing wave angular momentum, and can result in the planet spiralling 
into the star. However, in previous simulations the waves have not been observed to 
undergo instability for smaller amplitudes. Here we perform a stability analysis of a 
nonlinear standing internal gravity wave in the central regions of a solar-type star. 
This work has two aims: to determine any instabilities that set in for small- amplitude 
waves, and to further understand the breaking process for large- amplitude waves that 
overturn the stratification. Our results are compared with the stability of a plane in- 
ternal gravity wave in a uniform stratification, and with previous work by Kumar & 
Goodman on a similar problem to our own. Our main result is that the waves undergo 
parametric instabilities for any amplitude (in the absence of viscosity and thermal 
conduction). However, because the nonlinearity is spatially localised in the innermost 
wavelengths, the growth rates of these instabilities tend to be sufficiently small that 
they do not result in astrophysically important tidal dissipation. Indeed, we estimate 
that the modified tidal quality factors of the star that result are > 10^, and possi- 
bly much greater, which implies that the resulting dissipation is at least two orders of 
magnitude weaker than that which results from critical-layer absorption. These results 
support our explanation for the survival of all currently observed short-period planets 
around solar- type main- sequence stars: that planets unable to cause wave breaking at 
the centre of their host stars are likely to survive against tidal decay. This hypothesis 
will be tested by ongoing and future observations of transiting planets, such as WASP 
and Kepler. 

Key words: planetary systems - stars: rotation - binaries: close - hydrodynamics - 
waves - instabilities 



1 INTRODUCTION 

The tidal interaction between a short-period planet and its 
host star can result in evolution of the stellar and planetary 
spins and the planetary orbit. In particular, dissipation of 
the energy stored in the tidal response in the star can result 
in the planet spiralling into the star when the period of the 
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stellar spin is longer than the orbital period. This is because 
a final state in which the star spins synchronously with the 
planetary orbit cannot be achieved when the angular mo- 
mentum of the orbit is at most comparable with that of the 
stellar spin (Counselman 1973; Hut 1980). In addition, the 
star constantly loses spin angular momentum as a result of 
magnetic braking, which means that a synchronous equi- 
librium state does not exist, and the resulting tidal torque 
eventually acts to pull the planet towards the star (Barker & 
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Ogilvie 2009). The efficiency of tidal evolutionary processes 
depend on the dissipative properties of the star and planet, 
which are usually parametrised by a dimensionless quality 
factor^ Q' for each body, which is an inverse measure of the 
dissipation. This is usually defined to be proportional to the 
ratio of the maximum energy stored in a tidal oscillation to 
the energy dissipated over one cycle (e.g. Goldreich & Soter 
1966). The mechanisms that contribute to Q' for fluid bodies 
are poorly understood, but it is thought that Q' depends on 
tidal frequency, the internal structure of the body, and, in 
some cases, the amplitude of the tidal forcing. In this paper, 
we study the mechanisms of tidal dissipation in solar-type 
main-sequence stars, continuing an investigation described 
in previous work by the authors (Ogilvie & Lin 2007, here- 
after OL07; Barker & Ogilvie 2010, hereafter BOlO; Barker 
2011, hereafter Bll). 

The response of a fluid body to tidal forcing can be de- 
composed into a prolate spheroidal quasi-hydrostatic bulge, 
referred to as the equilibrium tide, and a residual wave-like 
response, which results from the nonzero forcing frequency 
in the frame of the fluid, often referred to as the dynami- 
cal tide. In radiation zones of solar-type stars, the dynam- 
ical tide takes the form of internal (inertia) gravity waves 
(IGWs) , which propagate at frequencies below the buoyancy 
frequency N. These have previously been proposed to con- 
tribute to Q' for solar-type stars (e.g. Goodman & Dickson 
1998, hereafter GD98; Terquem et al. 1998). 

A short-period planet efficiently excites IGWs at the 
top of the radiation zone (hereafter RZ), where their exists 
a location at which N ~ 1/P, with P being the planetary 
orbital period. These waves propagate downwards into the 
RZ, until they reach the centre of the star, where they are 
geometrically focused and can become nonlinear. If their am- 
plitudes are sufficiently large for the wave to overturn the 
isentropes^, the wave breaks and deposits its angular mo- 
mentum to form a critical layer, at which ingoing waves are 
efficiently absorbed. This results in a strong tidal torque, 
which can prevent the survival of sufficiently massive short- 
period planets around solar-type stars. However, it only oc- 
curs if the planet is sufficiently massive, or the centre of the 
star is sufficiently stably stratified. None of the planets cur- 
rently observed to orbit solar-type main-sequence stars is 
likely to excite waves that break, which could be an impor- 
tant explanation for their survival. 

In BO 10 and Bll we performed two- and three- 
dimensional simulations of these waves as they approach the 
centre of the star. The results were found to be very simi- 
lar in both two and three dimensions, and if the amplitude 
is insufficient for the waves to overturn the isentropes, the 
waves were observed to reflect approximately perfectly from 
the centre of the star, and no instability appeared to set in. 
However, this could be a result of insufficient spatial res- 
olution or run-time in the simulations performed thus far. 
In this paper we perform a detailed stability analysis of a 
standing internal gravity wave in two dimensions. The aim 

^ Related to the traditional Q by Q' = 3Q/2k, where k is the 
second-order potential Love number of the body. 

In stars, the stratification is actually composed of both entropy 
and composition gradients. When we refer to "isentropes" we ac- 
tually mean stratification surfaces, however, these are usually ap- 
proximately the same. 



of this work is to determine whether any instabilities are 
likely to occur in reality for small-amplitude waves, which 
are unable to overturn the isentropes. If an instability exists 
for these waves, and if this results in efficient tidal dissipa- 
tion, then this could have important consequences for the 
survival of short-period planets around solar-type stars. 

1.1 Stability analyses of IGWs 

Many stability analyses of a plane IGW in Cartesian ge- 
ometry with a uniform stratification have been performed 
(e.g. McEwan & Robinson 1975; Meid 1976; Drazin 1977; 
Klostermeyer 1982). These indicate that a monochromatic 
propagating plane IGW is always unstable to parametric in- 
stabilities, whatever its amplitude, in the absence of viscos- 
ity and thermal (or compositional) diff'usion. In that prob- 
lem such analyses were made possible for flnite-amplitude 
(in addition to infinitesimal amplitude) waves because the 
solution is exact. This is a consequence of the fact that the 
wavevector k and velocity u satisfies k-u = 0, implying that 
the advective operator u • V annihilates any disturbance be- 
longing to the same plane wave. These stability analyses 
allow a detailed understanding of the initial stages of the 
breaking process for these waves (e.g. Drazin 1977; Kloster- 
meyer 1982; Lombard & Riley 1996; or the review: Staquet 
& Sommeria 2002). 

When a small perturbation is added to a basic plane 
wave, the resulting evolutionary equations have periodic co- 
efficients. This allows the possibility for parametric instabil- 
ity to occur. The first study of this problem was by McE- 
wan Sz Robinson (1975), who considered perturbations with 
length scales much smaller than the primary IGW wave- 
length, in which case the problem can be reduced to the 
solution of Mathieu's equation. The motion of the fluid in 
the basic wave gives rise to unstable modes, just as para- 
metric oscillations of a pendulum are excited by periodic 
changes of its length. The growth rates of these parametri- 
cally unstable modes increase (linearly) with the amplitude 
of the basic wave. 

Subsequent analyses expanded the perturbation onto 
a Floquet basis, and relaxed the small-scale assumption. 
These studies all found that, in a dissipationless fluid, the 
disturbances with the largest growth rates have the small- 
est spatial scales (e.g. Drazin 1977; Klostermeyer 1982). In 
viscous or radiative fluids, dissipative effects scale with the 
inverse square of the length scale of a given mode. This 
means that the most unstable wavelengths will no longer be 
those of the smallest spatial scale, but will be those for which 
the competing effects of dissipation and (nonlinear) growth 
favour the latter, and this will depend on the Reynolds num- 
ber (also the Prandtl number when thermal diffusion is in- 
cluded) . 

Lombard & Riley (1996) & Sonmor & Klaassen (1997) 
performed a detailed stability analysis of a plane IGW, both 
demonstrating that the instability that contributes to wave 
breaking is driven by a combination of wave shear and wave 
entropy gradients. They find that wave- wave resonance in- 
teractions are the primary mode of instability for small- 
amplitude waves, with the picture being much more compli- 
cated near overturning amplitudes. However, no difference 
in the source of free energy driving the instability is found 
for waves that do and do not overturn the stratification for 
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some wave phase. Some of their results have been confirmed 
in recent high resolution numerical simulations (e.g. Fritts 
et al. 2009). We discuss these stability analyses further in 
relation to our results in ^8.3. 



1.2 This work 

In our problem we have obtained an exact 2D standing wave 
solution in cylindrical geometry representing IGWs near the 
centre of a solar-type star. This enables us to perform a sta- 
bility analysis of this wave for any amplitude. This is the 
subject of the present paper. One important difference be- 
tween our problem and previous studies is that the non- 
linearity is spatially localised to the innermost wavelengths, 
whereas for the plane IGW problem, though the nonlinearity 
may be localised within each wavelength, there is periodic 
repitition in space. 

In the centre of a star, (molecular) viscous damping is 
negligible, and the dominant linear dissipation mechanism 
is radiative diffusion. However, the waves excited by planets 
orbiting solar-type stars with several-day periods have much 
larger frequencies than their radiative damping rates. This 
means that parametrically excited modes with scales shorter 
than the primary wave could be produced. These will be 
damped by diffusion themselves, but not before they can 
draw energy from the primary wave, and possibly contribute 
to wave breaking. 

We have already demonstrated through direct numeri- 
cal simulations (BOlO; Bll) that a wave with sufficient am- 
plitude to overturn the stratification undergoes a rapid in- 
stability (with a growth time on the order of a wave period) 
which leads to wave breaking. We found that the wave over- 
turns the stratification during part of its cycle if the angular 
velocity in the wave exceeds the angular pattern speed of 
the forcing, i.e. it^/r > cj/m, where uj is the wave frequency 
and m is its azimuthal wavenumber. This can be expressed 
as the following breaking criterion derived in Bll valid for 
the current Sun. The tidally excited waves break near the 
centre of the star if the dimensionless nonlinearity in the 
wave 



;0.28 
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(1) 



where m^/p is the mass of the star /planet, and C is de- 
fined such that N — Cr near the centre of the star. The 
parameter A is defined so that the wave overturns the isen- 
tropes for some location in the wave if A ^ 1 (this is defined 
consistently with Eq. 16 below). This means that a one- 
day Jupiter-mass planet is not likely to excite IGWs with 
sufficient amplitudes to cause breaking at the centre of the 
current Sun. However, there is a strong dependence on C, 
which measures the strength of the stratification near the 
centre of the star, which implies that breaking is more likely 
in older and more massive stars (of solar-type, with radiative 
cores) . 

In the 2D simulations, the wave reflects perfectly from 
the centre of the star if its amplitude is insufficient to satisfy 
the breaking criterion, and long-term integrations do not 
show that any instabilities act on the waves. The picture 
in 3D is very similar. In this paper we perform a weakly 
nonlinear stability analysis of our 2D wave solution using a 
Galerkin spectral method. This work has two main aims: 



• to better understand the early stages of the breaking 
process for large-amplitude waves {A > 1); 

• to determine what (if any) instabilities may set in for 
waves that are unable to overturn the isentropes at any lo- 
cation in the wave {A < 1). 

The motivation for this study is that if the waves are 
subject to parametric instabilities (as proposed by Good- 
man & Dickson 1998, hereafter GD98; Kumar & Goodman 
1996; hereafter KG96), whatever their amplitudes, the re- 
flection of waves from the centre of the star will not be per- 
fect. This would stand in contrast to the prediction from 
linear theory, and the results of our numerical simulations. 
The simulations performed thus far may not have the spatial 
resolution or have long enough run time to be able to cap- 
ture small-scale parametric instabilities. Alternatively, the 
adopted boundary conditions may exclude the existence of 
parametric instabilities. If they indeed occur in reality, and 
the tidally excited waves are weakly nonlinearly damped by 
parametric instabilities, this could contribute to the tidal 
dissipation, and have implications for the survival of short- 
period planets with insufficient masses to satisfy Eq. 1 and 
cause breaking. 

The paper is structured as follows. First, we present 
the Boussinesq-type model derived in BO 10, and obtain an 
exact wave solution that represents a standing IGW that 
is confined within a circular domain. We then derive the 
equations governing linear perturbations to this wave, and 
expand these perturbations using an appropriate, complete, 
set of basis functions. The resulting eigenvalue problem is 
solved both for cases in which the wave does and does not 
overturn the isentropes, and the properties of the resulting 
unstable modes are studied. We particularly concentrate on 
determining the growth rates of the instabilities and how 
these vary with the parameters of the problem, as well as un- 
derstanding what is the source of free energy driving them. 
This is then followed by a discussion, in which we compare 
our results with previous studies of the stability of a plane 
IGW. We also compare our results with previous work by 
KG96 on parametric instabilties of tidally excited waves, and 
discuss the implications of our results for the tidal dissipa- 
tion in solar-type stars, and to the survival of short-period 
planets in orbit around them. 



2 INTERNAL GRAVITY WAVE STABILITY 
ANALYSIS 

We start with the adiabatic Boussinesq-type system (BO 10) 

D\i = -Vq + rb, (2) 

Db + C^r'U = 0, (3) 

V • u = 0, (4) 

D = dt + u-V, (5) 

where u is the fluid velocity, 6 is a buoyancy variable (pro- 
portional to the entropy perturbation) and ^ is a modified 
pressure variable. These equations were derived in BO 10 
from the equations of gas dynamics and are able to describe 
the dynamics of nonlinear IGWs near the centre of a star 
where the density is nearly uniform and the buoyancy fre- 
quency is proportional to radius. In this model N = Cr, 
where C is a constant that measures the strength of the 
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stable stratification at the centre. This model is valid in 
the innermost < 3% of a solar-type star, which contains 
multiple IGW wavelengths for the waves excited by short- 
period planets. Acoustic waves have been filtered out from 
this model. We have also omitted viscosity and thermal con- 
duction from these equations, although we will add these 
effects later. 

Since we restrict our problem to two dimensions, we 
can express the velocity field in terms of a streamfunction 
defined in polar coordinates (r, 0) by 



Ur 



-drip, 



(6) 
(7) 



which automatically enforces the solenoidality constraint on 
the velocity. We consider a circular region with r G [O,rotit], 
taking rout = 1, which is an impermeable outer boundary at 
constant entropy, i.e., ip{l,(j),t) = 6(1, (/),t) = 0, to confine 
the modes. We also adopt an inner regularity condition at 
r = 0, which chooses the regular solutions of the system^. 
This choice of boundary conditions ensures that the total 
energy of the perturbations is conserved, since the energy 
flux through the boundaries is always zero (because Ur = b = 
at r = 0, 1). We use dimensionless units such that the unit 
of length [L] — rout, the unit of time [T] — N~^^ = C~^r~^^^ 
and hence C = 1 in these units. 

To eliminate the modified pressure variable ^, we take 
the curl of the momentum equation: 



dt{V X u) = V X r6 - V X (u • Vu). 



(8) 



The ^-component of this equation gives the vorticity equa- 
tion, which, together with the buoyancy equation, is 



+ d4.h 
dtb + 



J{^,b), 



(9) 
(10) 



where the vorticity is C = —V'^ip. The nonlinear terms have 
been written in the form of Jacobians, defined by 



J{A,B) 



1 d{A,B) 
r d{r, (j)) 



1 



[{drA){dc^B)-{dc^A){drB)]. (11) 



We consider a stationary, stably stratified background 
containing a nonlinear wave (denoted by subscript w)^ sub- 
ject to a perturbation (denoted by primes). That is, we ex- 
pand 



(12) 
(13) 



The linearisation of Eqs. 9 and 10 in terms of the perturba- 
tion is 

at(-VV) + ^0^' = J(V^^,-VV) + -^(V^',-VV-) (14) 
dtb' + d^iP' = J{il;^,b') + J{il;',bn,), (15) 

which is two equations for two unknowns (V^^, b'). The non- 
linearities in this system provide coupling between differ- 
ent waves. We neglect the terms J('0^ — V'^ip') and J('0^ b')^ 
which is consistent with our weakly nonlinear approach. 



^ However, in the computation of the table of integrals described 
below, this is replaced by an impermeable inner boundary at 
rin = 10 to avoid the coordinate singularity at the origin. 
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with an arbitrary amplitude for illustration. Contours of constant 
ipw are the streamlines of the primary wave flow. The flow goes 
clockwise around the red streamlines and anticlockwise around 
the blue streamlines. Stagnation points are located at the radial 
nodes r where J2(kr) = 0, at azimuthal locations ^ = (2?tH-i) ^ 
for n G Z. 



2.1 Exact primary wave solution in 2D 

We consider a nonlinear gravity wave (the primary wave) 
with a well-defined angular pattern speed and azimuthal 
wavenumber m — 2. Eqs. 9 and 10 are invariant under trans- 
formation to a rotating frame, if the streamfunction is trans- 
formed appropriately. This is because the Coriolis force can 
be written as the gradient of a potential, and therefore has 
no effect. In the frame in which the wave is steady and <j) is 
the azimuthal coordinate, our primary wave is 



Re \ —AJ2{kr)e^ 



2_ 

kipw 5 



(16) 
(17) 



which is an m = 2 wave with Up radial nodes (to be cho- 
sen later), where J2 is a Bessel function. In general, A G C, 
but is time-independent in this frame. From here on, we take 
A G M, without loss of generality. Note that V'^tpw = —k^ipw, 
which implies that this solution is an exact (nonlinear) so- 
lution of Eqs. 9 and 10. We choose k such that J2{k) — 0. 
This is equivalent to confining the primary wave in a circular 
region of unit radius with an impermeable outer boundary 
at constant entropy. We plot an example of this wave with 
rip = 2 in Fig. 1. 

This wave overturns the stratification when drS < 0, 
where s = {1 / 2)r^ -\- bw is proportional to the total entropy. 
This is equivalent to ^drbw < —1- Note that overturning oc- 
curs only when A > 1, and is more likely for waves with large 
radial node numbers and small azimuthal wavenumbers. The 
size of the convectively unstable region of the m — 2 primary 
wave can be illustrated for a given A and n^, by calculating 



gdrS = r{r + drb) 



2 2A 

r H — — r (Ji(/cr) — J3(/cr)) cos 20. 
k 



For illustration, we plot the 2D region that is convectively 
unstable for several A values when rip = 2 in Fig. 2. An 
approximate size for the overturning region for small r when 
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Galerkin basis. This basis is adopted for two reasons: the 
Unear solutions (A — 0) take the same form, and it auto- 
matically satisfies the chosen boundary conditions. 

Note that our spectral-space amplitudes 
'ipm,n{t),bm,n{t) G C, SO wc must take the real part at 
the end of the calculation to obtain physical quantities. 
For each m, there is an infinite number of components 
with different values of n. In our spectral representation 
of the solution, we truncate these infinite series such that 
1 — Lm ^ m ^ Lm — 1, where Lm is an odd number, and 
^ n ^ Ln. This truncation is chosen so that we have an 
exactly equal number either side of m = (and a similar 
number either side of the primary wave m — 2) which 
ensures that our mathematical realisation of the problem 
has the symmetry property that we discuss in § 5.3 below. 



m m - 



-2.5 



Figure 2. Spatial exent of the region that is made convectively 
unstable by the primary wave for A = 1.1 (top) and 10 (bottom), 
with Up = 2. This region expands from the point r = when 
A = 1 to encompass the innermost few wavelengths for larger A. 



A> lis 



(18) 



at the most unstable wave phase. When A — 1 the over- 
turning region is the point r = 0, with the region expanding 
for larger A. If the instabilities that cause wave breaking are 
convectively driven, we would expect them to be strongly 
(though not necessarily completely) localised within these 
convectively unstable regions of the primary wave. 



2.2 Infinitesimal perturbations 

We consider linear perturbations to this finite-amplitude pri- 
mary wave, which we expand as (dropping the primes from 
now on) 



oo oo 

= ^ ll^m,n{t)Jm{km,nr)e' 
m— — oo n—0 
oo oo 

b = ^ hm,n{t)Jm{km,nr)e''' 

m= — oo n=Q 



(19) 



(20) 



where /cm,n is chosen such that the solutions for each az- 
imuthal wavenumber m, and radial node number n ^ 0, 
satisfy the outer boundary condition at r = 1, for which 



Jm{km,n) — 0. 



(21) 



This condition forces km,n ^ IR, for |m|,n G Z^. The above 
expansion automatically enforces a regularity condition on 
the perturbations at r = 0. Eqs. 19 and 20 define our 



2.3 Derivation of the evolutionary equations 

Evolutionary equations for the amplitudes ipm,n{t) and 
bm,n{t) can be derived by projection through integration 
onto the Galerkin basis. To do this we substitute the above 
expansions into the linear system defined by Eqs. 14-15. An 
important orthogonality relation is 



Jo Jo 



tim — m 



00 



= 7T[Jm-\-l{km,n)f Sn,n'Sm,m' , (22) 



where S is the Kronecker delta. Note that this results in 
different normalisation factors for each m and n wave. Also 
note that 

V ^Jm{km,nT)e — kjj^j^Jm{km,nf'^G . (23) 



2.4 Linear solutions (in the absence of the 
primary wave) 

Consider Eqs. 14-15 with J(. ..,...) = 0, which is equiv- 
alent to having a hydrostatic background with no primary 
wave flow. If we substitute the expansions Eq. 19-20 into 
Eqs. 14-15, and then multiply by rJm'{km'n'T)e~'''^ and 
finally integrate over G [0, 2n] and r G [0, 1], we obtain 



+ imbm,n = 



(24) 
(25) 



for each m, n, after relabelling m' ^ m, and n ^ n follow- 
ing the integration. This system, together with the boundary 
conditions, can be solved to give 



(26) 



with Ayn^n,Bm,n G C, and cjm,n - ±m/km,n- TMs is the 
frequency of a non-interacting wave which is able to exist 
in the container in the absence of any primary wave flow. 
If we do not truncate the Galerkin basis at some finite val- 
ues of Ln and Lm, then these modes would be dense in the 
frequency interval (0, 1), since the maximum buoyancy fre- 
quency Nmax = 1. In a frame in which the fluid rotates with 
angular velocity Q, the Doppler-shifted wave frequency is 
Coyn^n — (^m,n — TTiQ. Notc that km,n iucrcascs with both m 
and n, but 0Jm,n decreases with n and increases with m. 
When substituting the above solution back into Eq. 19-20, 
we obtain the linear solutions of the system, which are Bessel 
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functions of a given order \m\ with n nodes in the radial di- 
rection. This motivated our choice of Galerkin basis. 



with the integrals 



2.5 Nonlinear terms 

We obtain our system of equations from Eqs. 14-15 through 
the same approach as in the previous section, to obtain for 
each m and n, 



+J{i;,-V^i;y,)] |c/r#, (27) 

27r pi 



bm,n + im^i 



7v[Jm+l{km,n)]'^ Jo Jo 

^)]} 

_ 1 r^" f 

7r[Jm-\-l{km,n)]'^ Jo Jo 

+J(^,6^)] |dr#, (28) 

where the Jacobians contain sums over n' and m' . The sum 
over m is reduced to a pair of terms through the (j) inte- 
gration, using Eq. 22. A set of coupling integrals of triple 
products of Bessel functions also results, for which there is 
a sum of such terms over n\ i.e., an m, n wave is coupled 
through nonlinear terms to waves with m dz 2 and (in princi- 
ple) all node numbers n G {0, . . . , oo}. The system reduces 
to 

°° r 

+ imbm,n = ^ S {km-2,n' — k2,np) ^i^m-2,n' 

n'=0 ^ 

) A*V^^+2,n' (29) 

oo f 

bm,n + im'lpm,n = X! ) ^^rn,n,n' ^ (^m-2,n' " ^2,np '0m-2,n' ) 
n'=0 ^ 

^* {hm+2,n' - k2,np'lpm+2,n') |, (30) 



where 



(31) 



The coupling coefficients are 
2 



^m,n,n' 



[Jm-\-l{km,n^] 
^ V 

normalisation 

2 

[Jm + l(k m,n I \ 
normalisation 



((m 2)X^^^^^/ 2X^^^^^/) , (32) 



((m + 2)I^,„,„, +2I^,„,„,), (33) 



-^m,n,n' 


/ Jm{km. 

Jo 


,nr) [drJ2(k2,npr)] Jm-2{km-2,n' r)dr, 


, - 


/ Jm{km. 

Jo 


,nr)J2{k2,npr) [drJm-2{km-2,n'r)] dv, 


, = 

-^m,n,n' 


1 Jmikm. 

Jo 


,nr) [drJ2(k2,npr)] Jm+2{km+2,n' r)dr, 


, - 


/ Jmikm. 

Jo 


,nr)J2{k2,npr) [drJm+2{km+2,n''r)] dr. 



Note that these are related by 

^m,n,n' ^m-\-2,n' ,n •> (34) 

■^m,n,n' ■^m-\-2,n' ,n ■^m-\-2,n' ,n 5 (35) 

where the latter follows from an integration by parts. For 
use in the derivation of the spectral space energy equation 
in a subsequent section, we find it convenient to define 

^m,n,n' — TT [ Jrn,+1 (/Cm,,n)] OLm,n,n' ^ (36) 

and similarly for Prn,n,n'- This is because we then have the 
relation 



^m,n,n' ^m+2,n',n' 

2.6 Diffusive terms 



(37) 



In the presence of viscosity and radiative diffusion (or hy- 
perdiffusion) Eqs. 14 and 15 have the additional terms 
(-l)2+"zyV^+^"^ and (-l)^+^/^V^"6, respectively. Here a 
is chosen to give the standard diffusive operator {a = 1), 
or hyperdiffusion (a = 2, 3). In this case we obtain the (lin- 
earised) system 

+ imbm,n = -Z^A^mTn (38) 
^m,n ~h irniprn,!!, — l^k^^^bm,n') (39) 

instead of Eqs. 24-25, where u is the kinematic (hyper-) vis- 
cosity and K, is the thermal (hyper-) diffusivity, with similar 
modifications to Eqs. 29-30. The dispersion relation is then 

^2 



I ^m,n + i^kr^^ri) {^m 



h.2 



(40) 



indicating that the frequencies of the allowed solutions are 
modified in the presence of z/, The growth rate expected 
in the presence of weak diffusion is therefore Im [uj] — |(i/ + 
i^)k^^n^ with Im [cj] being the appropriate inviscid growth 
rate. Since Eqs. 29-30 allow instability, to obtain growing 
modes in the presence of diffusion, v and hi must be suf- 
ficiently weak so that diffusive terms do not dominate over 
the nonlinear terms except for values of n and m close to the 
resolution limits of Ln and Lm- Hyperdiffusion with a = 3 
is adopted since it better restricts the dissipation to the 
highest wavenumbers. This enables numerical convergence 
in the eigenvalue problem discussed below, but does not sig- 
nificantly perturb the growth rates of the lower wavenumber 
eigenmodes with the values of i/, n that we adopt. From here 
on, we also take u = k,. 

With the inclusion of diffusive terms we require 1 -\-2a 
additional boundary conditions at each boundary. These are 
regularity conditions at the centre, and at the outer bound- 
ary we can consider a variety of idealised boundary condi- 
tions of the form V^^^ = for a = 0, 1, . . . , a and V^^6 = 
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for a — 0, — 1. These are automatically satisfied 

by our Galerkin basis Eq. 19-20 and the definition of km,n- 
With these boundary conditions (which force km,n to be 
real) , the Galerkin basis is therefore exact for the single- wave 
diffusive problem. Note that this is also true if hyperdiffu- 
sion is adopted, due to property Eq. 23. In the numerical 
solution of the eigenvalue problem in the following sections 
we include these diffusive terms to achieve numerical conver- 
gence. This is necessary because in the absence of diffusion, 
the most unstable modes are found to prefer the smallest 
spatial scales. 



3 METHOD OF SOLUTION 

Our system Eqs. 27-28 can be written in the form of a gen- 
eralised eigenvalue problem of the form 



AU = cjBU, 



(41) 



where U is the column vector whose components are the 
quantities (-?/;, 



m,n 5 ^m,n ^ 



for each m and n. A is the block 
tridiagonal matrix representing the system. This is done by 
seeking normal mode solutions of the form ipm,n{t) oc e~*^* 
for each m and n, and similarly for bm,n{t). B is the diagonal 
matrix that can be composed as blocks of the form 







(42) 



for each m and n. We solve this problem using standard 
generalised eigenvalue solver routines, such as ZGGEV in 
the LAPACK library. This returns the eigenvalues {c<;}, and 
the spectral space eigenfunctions {ipm,n,bm,n} correspond- 
ing to each eigenvalue. The real space eigenfunctions can be 
reconstructed from these, using Eqs. 19-20. 

We choose Ln = 50 and Lm = 27 for most of the cal- 
culations (though a small number of higher resolution cal- 
culations were performed with Lm values up to 45, which 
confirm that our results are not dependent on resolution). 
There is a limit to the maximum values of Ln and Lm that 
we can reasonably adopt, due to the computational cost of 
choosing large values for each of these parameters. One rea- 
son for this is that the nonlinear terms require the compu- 
tation of a large number of integrals, many of which have 
highly oscillatory integrands, and require very small relative 
error tolerances to be computed accurately (by the method 
of computation that we will describe in the next subsec- 
tion). Another reason is that we compute the eigenvalues 
and eigenvectors using a QZ alogorithm, which has a high 
computational cost 0(105'^) for an 5* x 5* matrix, where 
S = 2Lm {Ln — 1). This limits the values of Ln and Lm- 
With our choice of a = 3 hyperdiffusion, it has been found 
that 10~^^ ^ ^ 10~^^ is appropriate. This hyperdiffusion 
is found to give the spurious eigenvalues, whose eigenfunc- 
tions oscillate at the smallest scales, and are therefore not 
adequately resolved, a large decay rate, and allows our grow- 
ing modes to be adequately converged for the values of A 
and Up that we consider. 



3.1 Numerical computation of table of integrals 

The tables of integrals defined in § 2.5 are computed for 
each value of m, n, using a 4th/5th order adaptive step 



Runge-Kutta integrator. To enable efficient computations, 
the Bessel functions are computed simultaneously with the 
integrals. Note that Bessel's equation 

drirdr^) + r - ^ j V = 0, (43) 

can be rewritten as the coupled set of first order ODEs 



^2 '^m,n I 



dr 

dr r 

We also need the derivatives of various Bessel functions, so 
we also integrate 



(44) 
(45) 



dr2 



(m^ - r^fc^,„)V' - (. 



(46) 



to obtain the first derivative of each Bessel function. To com- 
pute the integrals, we integrate the integrands of X^^^^^/ for 
i G {1,2,3,4}, and take the value at r = 1. This involves 
solving a system of 15 ODEs in total for each m,n,n^ Us- 
ing our chosen resolution of Ln = 50, Lm = 27, this involves 
the computation of Lm{Ln — 1)^ ~ 10^ integrals in total. 
These are computed for a given number of radial nodes in 
the primary wave in the range ^ rip ^ 12. We impose an 
inner boundary at r^n = 10~^ to avoid the coordinate sin- 
gularity at the origin, and use initial conditions appropri- 
ate from considering the asymptotic behaviour of the Bessel 
functions. For small Ln and Lm values the numerical inte- 
grals have been checked to agree with those computed from 
Mathematica, and for large Ln and Lm, several integrals 
containing the highest n and m values were also checked. 
We use a relative error tolerance of 10~^^, which has been 
found to compute the most oscillatory integrals (correspond- 
ing to the highest n and m value) accurately (compared with 
Mathematica) to within at least 6 decimal places. 



4 KINETIC AND POTENTIAL ENERGY 
EQUATIONS 

The kinetic and potential energies can be computed from 
either the real-space or spectral-space eigenfunctions. This 
enables a determination of the dominant source of free en- 
ergy driving the instability (e.g. Lombard & Riley 1996), and 
also provides an independent calculation of the growth rate, 
which can be used to check our numerical code. We derive an 
energy equation in spectral space, and compute the volume- 
integrated terms using the numerically computed eigenfunc- 
tions, without converting to real space. This has been found 
to reduce numerical errors, resulting from large numerical 
cancellations in the most oscillatory Bessel functions, when 
the energy equations are instead computed in real space. 
In addition, it is simpler to construct the (hyper-) diffusion 
terms in spectral space, so that they can be fully taken into 
account in the energy budget. 
We define 



'rdrd4>, 



b rdrd4>, 



E = K + P, 



(47) 

(48) 
(49) 
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as the kinetic, potential and total energy densities of the 
disturbance, respectively. Note that 



^ oo oo 

^ = 2 ^ Xl^m,n|'?/^m,n|^7r[Jm+l(A:m,n)]^ , (50) 
m= — oo n=0 
^ oo oo 

P = 2 X! ^^\^rn,n\'^7r[Jm+l(km,n)f , (51) 
m= — oo n=0 

E = K + P. (52) 

Evolutionary equations for the volume-integrated energy 
can be obtained from Eq. 29 and 30 together with the 
(hyper-) diffusion terms. After some rearrangement, these 
can be written 



E = J\fsw + Mhw + F^ + F^, 
where 



(53) 
(54) 
(55) 



OO oo oo 



m= — oo n=0 n'=0 



oo oo oo 



J\fbw = Re ^ ^ ^ iQ;^,n,n'^^2,np 
m= — oo n=0 n'=0 



oo oo 



= Re ^ ^ ^mV^m,n^m,n7r [ Jm + 1 (A^ri 
m= — oo n=0 
oo oo 

= Re ^ ^ -i^^^n"|^m,n|^7r [ Jm+1 (/i^m,n)] ^ , 

m= — oo n=0 

oo oo 

F^ = Re ^2 ^2 -^^^n|^m,n|^7r [Jm+l{km,n)f • 
m= — oo n=0 

JVsw represents the production of perturbation kinetic en- 
ergy from the primary wave shear. N'bw represents the pro- 
duction of perturbation potential energy from the primary 
wave entropy gradients. Whichever of JVsw or JVbw is domi- 
nant tells us whether this instability is driven by the wave 
shear or the wave entropy gradients. Fb is the buoyancy flux 
term, representing conversion between kinetic and potential 
energies of the disturbance. Finally, Fi^ and F^ represent the 
irreversible loss of kinetic and potential energies as a result 
of the (hyper-) diffusion. 

After truncation at I and each 

of these terms are computed from the spectral space eigen- 
functions, together with the numerically computed table of 
integrals. The growth rate can then be computed from 



Im[ct;'] = 



E _ K _ P 
2E ~ 2K ~ 2P' 



(56) 



We have checked that each of these equations are accurately 
satisfied to within at worst a few percent for each of the un- 
stable modes discussed in this paper. This provides a check 
of our analytical derivations and numerical calculations, and 
should convince ourselves that our results are consistent. 



5 NUMERICAL TESTS 

In this section we briefly mention several numerical tests 
which we have performed to validate our numerical code, 
in addition to the one mentioned in the previous section. 
Following this section, in § 6 and 7 we discuss the results 
of our stability analysis for waves with A < 1 and A > 1, 
respectively. 



5.1 Linear 

In the absence of nonlinear couplings (A = 0), we ob- 
tain a set of non-interacting modes with eigenfreqencies 
00 G {uJm,n}^ where uJm,n — -^m/km,n (in the inertial frame), 
as we predicted in § 2.4. In the absence of diffusion, these 
have zero growth rate, i.e., Im[cL;] = 0, for all eigenmodes. 
When hyper diffusion is included, the eigenmodes each have 
a nonzero decay rate determined by the values of v and At, 
which is very accurately (to more than 10 decimal places) 
computed from considering only the terms F^^ and F^ in 
the energy equation. The real-space eigenfunctions that re- 
sult are what is predicted from linear theory, in that they 
are Bessel functions of order m with n nodes in the radial 
direction. 



5.2 Weakly nonlinear 

We followed a typical eigenmode as A is gradually increased 
from zero, and found that for \A\ <C 1, the shift in the eigen- 
frequency Re[^C(;] oc A^, as we would expect for modes not 
undergoing parametric resonance (e.g. Landau & Lifshitz 
1969). 



5.3 Symmetries 

The real-space solutions can be represented in the form 



Re 



i{m4> — u!t) 



(57) 



This is symmetric under the transformations Cm,n{r) 
Cm.ni'^)^ —'^ £^nd uj — cj*. This symmetry should 

exist for all primary wave amplitudes, and results from the 
fact that only the sign of the pattern speed uo /m has mean- 
ing, and not the sign of the wavenumber or frequency. This 
means that when the eigenvalues are plotted on the complex 
frequency plane, they should be symmetric about Re [uS\ = 0. 



6 RESULTS FOR WAVES WITH A < 1: 
PARAMETRIC INSTABILITIES 

We examine the unstable modes that exist when < A < 1, 
which is when the primary wave does not overturn the strat- 
ification at any location in the wave. The instability is a 
parametric instability, for which a simple model is briefly 
reviewed in Appendix A. When A / 0, the fraction of eigen- 
modes that are growing is nonzero (above a critical A set 
by the values of Up and zy, which can be understood from 
Eq. A7) and increases with A, as nonlinear growth starts to 
dominate over the decay due to diffusion for a larger number 
of modes. For small A, the eigenvalues of the unstable modes 
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Figure 3. Distribution of unstable eigenvalues on the complex 
frequency plane for A = 0.1, Up = 2 and ly = 10'^'^, 10'^^, IQ-^^. 
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Figure 4. Distribution of unstable eigenvalues on the complex 
frequency plane for various values of A, with Up = 2 and u = 

10-14. 



displayed on the complex plane are distributed in two curves 
which are symmetric about Re [uj] =0. This is illustrated in 
Fig. 3, and is a result of the symmetry described in § 5.3. 

In the limit that A ^ 0, the unstable modes can be 
identified as the parametrically excited free- wave modes. 
These are a pair of free wave modes that exist when A = 0, 
which undergo modifications to their complex frequencies 
at 0{A) that reinforce each other. We have verified that the 
modes consist of a pair whose frequencies approximately add 
up to Up in the inertial frame, with a detuning /S^/uOp < 10~^. 
As A is increased, the unstable modes consist of gradually 
more complicated superpositions of free wave modes, un- 
til for A > 1, the eigenfunctions become localised in the 
convectively unstable regions. This will be studied in more 
detail in § 7. For A < 1, the eigenfunctions exist because of 
their confinement by the boundaries, though they interact 
quite strongly with the primary wave, and are generally not 
simply free wave modes with a nonzero growth rate. 

The number of unstable modes that exist in this am- 
plitude range depends quite strongly on viscosity. This is 
illustrated in Fig. 3. The number is also found to decrease 
as Up is increased. These two behaviours are related by 
the fact the a wave should undergo parametric instabili- 
ties, which have largest growth rates when the resonant tun- 
ing is good, which is more likely to occur for perturbations 
with larger wavenumbers. However, these large wavenumber 




-2.5 



-0.8 



-0.6 -0.4 

logio ^ 



-0.2 



Figure 5. Im[a;] vs. A for the most unstable mode with Up = 
2 and u = 10-^^,10-^^. The solid red line has slope 1. This 
shows that the instability approximately scales linearly with A 
for small A, indicating that the instability when A < 1 is due to 
a parametric resonance. 



components are strongly damped by diffusion. Increasing Up 
means that the "effective resolution" available to capture the 
unstable modes decreases. This is the same as increasing z/, 
hence the same trends exhibited in increasing Up and ly. 

The neat distribution of eigenvalues into two curves 
does not persist as A is increased, as illustrated in Fig. 4. 
The frequencies (in the rotating frame) are primarily smaller 
than the primary wave frequency in the inertial frame. How- 
ever, a small number of modes exist for A > 0.3 which have 
frequencies larger than ujp. Nevertheless, in each case the 
frequencies of the unstable modes are always smaller than 
the maximum buoyancy frequency in the fiow (which corre- 
sponds with l/cjp ^ 6 in the units of this figure, for rip = 2). 
This makes sense if these are gravity wave- like disturbances, 
which are parametrically excited by the primary wave. 

In Fig. 5, the growth rate for the most unstable mode 
is shown to scale approximately linearly with A for A ^ 1. 
A slope of 1 in this figure is predicted for A ^ 1 if the 
instability is due to a parametric resonance, and this appears 
to approximately hold for all A in this range. 

Our most important result of this study is illustrated in 
Fig. 6. This shows that the growth rate scales inversely with 
the number of wavelengths within the domain (note that 
this is after normalising by ujp). In this figure we plot the 
logarithm of the growth rate versus log^Q Up for A = 0.5 and 
0.8, respectively. The slope is always approximately equal to 
— 1 for low Up. The tail-off at larger Up is due to diffusion, 
and arises because modes with smaller spatial scales para- 
metrically excite modes with even smaller spatial scales (as 
a result of the theorem proved by Hasselmann 1967), which 
are then more easily damped by diffusion. As we would ex- 
pect from this interpretation, the value of Up at which dif- 
fusion dominates moves to smaller rip as ly is increased. The 
inverse dependence on Up that is present when diffusion is 
unimportant is a key result. This suggests that although 
parametric instabilities exist for any amplitude in the ab- 
sence of diffusion, in a sufficiently large domain they become 
unimportant. We discuss the relevance of this result to tidal 
dissipation in § 8. 
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Figure 6. logj^Qlmfcj/tJp] vs. log^o for the most unstable mode 
when A = 0.5 and 0.8 respectively, for u = 10"^^ and 10~^^. The 
red line has a slope —1. The tail-off at larger Up is due to diffusion. 



6.1 Eigenfunctions 

In the top panels of Figs. 7 and 8 we plot the real (blue solid 
lines) and imaginary parts (red dashed lines) along = 
of the spatial eigenfunctions of the most unstable mode for 
A = 0.1 and 0.5, respectively. We have taken rip = 2 and 
u = 10~^^. These modes exist throughout the box and are 
confined by the boundary. They become more distorted from 
the free wave modes as A is increased towards unity (and 
in fact also for A > 1, apart from the localised modes), es- 
pecially in the innermost wavelengths of the primary wave. 
Note that the amplitude of the eigenfunctions is arbitrary 
since we are solving an eigenvalue problem. The bottom pan- 
els of Figs. 7 and 8 show the spectral space eigenfunctions 
of the same unstable modes. We have normalised ipm,n and 
bm,n to their maximum values and taken the base 10 loga- 
rithm of each component to produce the figures. These show 
that growing modes for the chosen values of A are not simply 
a pair of free wave modes which are excited by the primary 
wave. They contain many n and m values localised around 
a particular n and m, and are therefore interacting strongly 
with the primary wave. Multiple n and m values are involved 
even when A = 0.1. For the value of u adopted, these modes 
are well resolved, as is shown from the amplitude decay of 
the spectral space eigenfunction, which occurs before the 
resolution limit is reached in n and m. 



6.2 Energetics of the instabilities 

When A < 1, the isentropes are never overturned by the 
primary wave, so a pure radial convective instability is not 
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Figure 7. Top: Real (blue solid lines) and imaginary (red dashed 
lines) parts along </) = of the spatial eigenfunctions for the 
most unstable mode for A = 0.1, Up = 2, u = 10"-*^^. The 
eigenfrequency is uj/up = 0.378 + O.OOli. Bottom: Spectral space 
eigenfunction of the same mode. The colour scale represents 
logio |^m,n/max{'0rn,n}|, and similarly for bm,n- 



possible. However, these instabilities could be driven by the 
free energy resulting from the primary wave shear or entropy 
gradients, or a combination of the two. In this section, we 
compute the spectral-space energy contributions outlined in 
§ 4, for a sample of growing modes in this amplitude range. 
We have confirmed that the growth rate is accurately com- 
puted from Eq. 56, to within a few percent at most. In Ta- 
ble 1 we outline the contributions to the growth rate from 
each term in Eqs. 53-55 for the most unstable mode for 
A = 0.1,0.5 and 1 with = 2 and ly = 10"^^. The eigen- 
functions for two of these are plotted in Figs. 7 and 8. These 
examples are illustrative of every unstable mode that exists 
when A < 1 (and also the non-localised modes that exist 
when A> 1). 

Firstly, we note that the integrated kinetic and poten- 
tial energy of the modes are in approximate equipartition. 
A single wave can be proved to be in exact equipartition 
(as is shown in Appendix B), so we would expect K ^ P if 
these modes are parametrically excited gravity waves with a 
single n and m. That they are in approximate equipartition 
and include many n and m components indicates that these 
modes are the larger A generalisations of the parametrically 
excited free wave modes. The source of free energy driving 
these modes is entirely the potential energy resulting from 
primary wave entropy gradients. Somewhat surprisingly, the 
primary wave shear contribution is much smaller, and actu- 
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Figure 8. Top: Real (blue solid lines) and imaginary (red dashed 
lines) parts along </) = of the spatial eigenfunctions for the 
most unstable mode with A = 1, rip = 2, u = lO"-*^^. The 
eigenfrequency is ou/oup = 0.540 + O.OlTz. Bottom: Spectral space 
eigenfunction of the same mode. The colour scale represents 
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Figure 9. Distribution of unstable eigenvalues on the complex 
frequency plane for Up = 2 and various A with i/ = lO"-*^^. 
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Figure 10. Distribution of unstable eigenvalues on the complex 
frequency plane for A = 10 and various rip, with u = 10"-*^^. 





A = 0.1 




A = 0.5 




A = 1 


Re [uj] 


0.065 




0.080 




0.093 


Im [uj] 


1.86 X 10- 


4 


1.97 X 10- 


3 


3.00 X 10-3 


K 


3.28 X 10- 


2 


3.82 X 10- 


2 


3.70 X 10-2 


P 


3.25 X 10- 


2 


3.87 X 10- 


2 


4.16 X 10-2 




-3.20 X 10 


-7 


-3.98 X 10 


-5 


-9.55 X 10-5 




3.57 X 10- 


5 


3.76 X 10- 


4 


6.62 X lO-'^ 


Fb 


1.81 X 10- 


5 


2.04 X 10- 


4 


3.51 X 10-4 




-5.58 X 10 


-6 


-1.35 X 10 


-5 


-3.29 X 10-5 


F^ 


-5.55 X 10 


-6 


-2.01 X 10 


-5 


-6.03 X 10-5 



Table 1. Energy components of the most unstable mode for A = 
0.1, 0.5 and 1, with Up = 2 and u = IQ-^^ . Note that oup ^ 0.17. 



ally stabilises the modes. This instability converts primary 
wave potential energy to disturbance potential energy, and 
then converts approximately half of this input energy to the 
disturbance kinetic energy through the buoyancy flux term. 
This process results in approximate equipartition between 
K and P. Note that the entropy gradients in the primary 
wave are insufficient to cause convective instability. These 
modes are driven by weaker entropy gradients in radius and 
azimuth. 



7 RESULTS FOR WAVES WITH A > 1: THE 
INITIAL STAGES OF WAVE BREAKING 

When A > 1, the primary wave overturns the stratification 
during part of its cycle. Our simulations in BO 10 have shown 
that an instability breaks the wave within a few wave peri- 
ods once this first occurs. The initial stages of this breaking 
process are examined in this section by choosing A > 1. To 
resolve convectively unstable modes with the adopted values 
of Ln and Lm , we require the size of the overturning region 
to be sufficiently large. Since overturning occurs only at the 
point r = when A = 1, this necessitates choosing values 
of A larger than unity to capture such instabilities. We are 
interested in instabilities which act to break the waves in an 
(effectively) unbounded domain (the central regions of the 
RZ of a solar-type star), therefore the appropriate unstable 
mode should not rely on the boundaries for confinement, 
and should be localised within the innermost wavelength of 
the primary wave. This is because the presence of confining 
boundaries is artificial, and is imposed to specify the prob- 
lem. With this in mind, we now discuss the results of our 
stability analysis for waves with A > 1. 

The eigenvalues of the unstable modes displayed on the 
complex plane are shown in Fig. 9 for A = 5 and A = 8, 
both with rip = 2 and v — lO"^"^. The most unstable modes 
are located on distinct branches, which stand above the 
continuation of the modes that exist when A < 1. From 
studying their eigenfunctions, we find that the modes on 
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Figure 11. lm[uS\ vs. A for the most unstable mode for rip = 2, 
i/ = 10-13 and 10-14. The instabihty grows within a primary 
wave period when A > 5, which is when the locahsed modes 
begin to appear. 

the branches are localised disturbances, unlike those below 
the main branches. The modes on the branches could there- 
fore represent the type of mode that breaks the primary 
wave (this is discussed further in the next subsection) . These 
branches extend further from the origin the greater the value 
of A. In Fig. 10 we plot the unstable modes for A = 5 for 
two values of Up. This shows that the unstable modes on 
the branches do not move around significantly as Up is var- 
ied. They therefore depend only weakly on the location of 
the outer boundary. Note, however, that the growth rate be- 
comes smaller as we go to larger Up because of the increasing 
importance of diffusion. 

Note that the largest frequency of some growing modes 
is larger than the maximum buoyancy frequency Nmax — 1 
(which corresponds with l/cjp ^ 6 in the units of this figure, 
for rip — 2). The nonzero frequencies of the modes in this 
frame indicate that they are oscillatory, and are non-steady. 
In addition, the growth rates of the most unstable modes are 
sufficiently fast compared with the primary wave frequency 
that the instability grows within several wave periods after 
onset. This is consistent with the results of our simulations 
described in BOIO. 

The growth rate of the most unstable modes for a given 
rip increases with A as illustrated in Fig. 11 for rip = 2, 
where curves for v — 10~^^ and 10~^^ have been plotted. 
There is an approximate square root dependence for A > 4. 
If the instability is driven by convectively unstable entropy 
gradients, then we might expect 

Im[cj'] < Vmax [-N'^] 

^maxj - - ^?^r (Ji(A:2,npr) 

>|X 1/2 
-J3{k2,npr)) COS 20 > j 

Thus, for large A the growth rate should scale with the 
square root of the primary wave amplitude. This behaviour 
is not observed when 1 < A < 4. In this range, the square 
root dependence may not be exhibited partly because there 
is insufficient resolution to accurately capture the modes 
that contribute to breaking since the overturning region is 
small compared with the box size. We have noticed that the 
growth rate does not significantly depend on u (and there- 




logio 

Figure 12. Im[a;] vs. rip for the most unstable localised modes 
when A = 5 and u = 10- and 10- 

fore the resolution), for the most unstable modes, except in 
the range 1 < A < 4, which supports this explanation. 

The behaviour of the growth rate on rip is illustrated 
in Fig. 12 for the most unstable mode when A — 5, for two 
values of v. For large rip, the unstable modes have sufficiently 
small spatial scales for diffusion to become important, so we 
expect a tail-off at large rip. The important point that can 
be taken from this figure is that the (normalised) growth 
rate of the localised modes on the branches does not depend 
on the number of wavelengths within the domain, for modes 
that are not strongly affected by diffusion. This means that 
the instability can be important in a large domain, such as 
a the RZ of a solar- type star, which contains many primary 
wavelengths. Note that this is very different to the behaviour 
found for the excited modes when A < 1, as shown in Fig 6. 

7.1 Eigenfunctions 

In the top panel of Fig. 13 we plot the real (blue solid lines) 
and imaginary parts (red dashed lines) along = of the 
spatial eigenfunctions of the most unstable mode for A = 5, 
rip — 2 and f — 10~^^. The spectral-space eigenfunction 
of this mode is plotted in the bottom panel. The mode is 
strongly nonlinear ly interacting with the primary wave, as 
shown by the number of different n and m values that appre- 
ciably contribute. The contribution to the eigenfunction is 
nonzero, but not maximal, near \nn\ — Lm 

1 , and is negligi- 
ble at n = Ln, which indicates that this mode is adequately 
resolved. The eigenfunction is spatially localised within the 
innermost wavelengths of the primary wave. Each of the sev- 
eral most unstable modes in the range 5 ^ A ^ 10 which 
lie on the branches in Fig. 9 are localised modes, and have 
qualitatively different form to the type of modes that exist 
below the branches, which are a continuation of the modes 
that exist when A < 1. As we go to larger A for the same 
value of z/, the most unstable mode utilises an increasing 
number of n and m values up to the resolution limit. This 
means that to adequately resolve the modes we would have 
to either increase the resolution or the value of v. 

The components of the spatial eigenfunction of the most 
unstable mode when A = 5, = 2 and v — 10"^"^, is plotted 
in Fig. 14 on the (0, r)-plane, to further illustrate the spatial 
dependence of this mode. In Fig. 15 we plot the region of 
negative A/"^ for the same primary wave. A comparison of 
these figures makes clear that the eigenfunction is primarily 
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Figure 13. Top: Real (blue solid lines) and imaginary (red 
dashed lines) parts along </) = of the spatial eigenfunctions 
for the most unstable mode with A = 5, rip = 2, = lO"-*^^. 
The eigenfrequency is uo/uop = 9.68 + 0.302i. Bottom: Spectral 
space eigenfunction of the same mode. The colour scale repre- 
sents logio |'0m,n/max{'0rn,n}|, and similarly for 6m, n- 



localised within the regions made convectively unstable by 
the primary wave entropy perturbation. This adds further 
evidence to the conjecture that the instability is convective. 
We also find that any unstable mode on the branches of 
Fig. 9 that are excited when 5 ^ A ^ 10 are similarly lo- 
calised and have a qualitatively similar appearance to the 
eigenfunction plotted in Fig. 14. 



Figure 14. 2D Spatial eigenfunction of the most unstable mode 
for A = 10, rip = 2 and u = lO"-*^"^, plotted on the ((/>, r)-plane. 
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7.2 Energetics of the instabilities 

In this section, we compute the spectral space energy con- 
tributions outlined in § 4 for a representative sample of the 
localised growing modes that exist when A > 1. We have 
confirmed that the growth rate is accurately computed from 
Eq. 56, to within a few percent for the modes considered 
in this analysis. However, it must be noted that the most 
unstable mode when A > 5 is typically not fully resolved 
with our adopted resolution and zy, in that there is nonzero 
power in the highest n and m values. This can lead to er- 
rors in the energy analysis typically of order 10 — 30%, so we 
leave these modes out of this analysis, and only choose those 
that are adequately resolved^ for Table 2. In this table, we 



Figure 15. Unstable region for A = 10, Up = 2. This can be 
compared with Fig. 14. 



outline the contributions to the growth rate from each term 
in Eqs. 53-55 for several unstable modes that exist when 
A > 5, each with rip — 2 and f — 10"^*^. The eigenfunction 
corresponding to the first of these is plotted in Figs. 13 and 
14. 

As in the case of the modes that exist when A < 1, the 
instability is driven by the free energy associated with pri- 
mary wave entropy gradients, as is shown by the fact that 
J\fhw is the dominant contribution to the growth. This is 
indeed what would be expected of a convectively driven in- 



^ A small number of higher resolution calculations with Lm val- 
ues up to 45 have been performed to fully resolve these modes. 
These calculations have confirmed that although the numerical 



values for several quantities may differ slighty for the most poorly 
resolved modes, our main results are not dependent on resolution. 
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A = b 


A = 7 


Re [uj] 


1.67 


2.40 


Im [uj] 


0.052 


0.134 


K 


6.95 X 10-2 


3.37 X 10-2 


P 


1.73 X 10-2 


1.02 X 10-2 


sw 


-5.50 X 10-3 


-1.57 X 10-3 




2.05 X 10-2 


1.55 X 10-2 


Fb 


1.82 X 10-2 


1.22 X 10-2 




-2.10 X 10-3 


-1.33 X 10-3 


F^ 


-0.517 X 10-3 


-0.80 X 10-3 



Table 2. Energy components of the most unstable mode for A = 
5 and 7, with Up = 2 and u = 10-^3 ^ote that ujp ^ 0.17. 

stability. In addition, the primary wave shear is much weaker 
and tends to stabilise the modes. Unlike the modes that ex- 
ist when A < 1, we do not necessarily have K ^ and 
examples have been found that do and do not satisfy ap- 
proximate equipartition, so these modes do not all appear 
to be gravity wave-like, unlike the parametrically excited 
modes that exist when A < 1. 

The growth rates are always < Re [V— , which is ex- 
pected to be an upper limit if the instability is convective. 
The negative contribution of shear, as well as hyper diffusion 
mean that the modes that we have calculated have some- 
what smaller growth rates that this simple estimate would 
predict. The route of energy transfer that drives the instabil- 
ity is the same as for the parametric instabilities discussed 
in the previous section. 



8 SUMMARY AND DISCUSSION OF RESULTS 

In the previous two sections we have analysed the instabili- 
ties that exist when A < 1 and A > 1. 

8.1 Wave breaking 

When A > 1 we have identified a class of localised modes 
that are driven by convectively unstable entropy gradients 
in the primary wave. These modes exist in the absence of 
an outer boundary, and are very likely to have initialised 
the wave breaking process in the simulations presented in 
BOlO and Bll. Subsequent stages in the breaking process 
are not studied using this stability analysis because these 
would involve nonlinear interactions between the perturba- 
tions to the wave, which we neglected in our weakly nonlin- 
ear approach. The instability growth time is of the order of 
a primary wave period, which is in agreement with the wave 
breaking times observed in our simulations. 

8.2 Parametric instabilities 

When A < 1, there exist pairs of parametrically excited 
modes driven by (convectively stable) primary wave entropy 
gradients, with wave shear playing a subordinate stabilis- 
ing role. These modes exist because of their confinement 
by the outer boundary. Our most important result regard- 
ing these modes is the inverse dependence of the growth 
rate on rip. This can be explained by considering the rela- 
tive time the primary wave spends in the innermost regions. 



where its nonlinearity is strongest. The fraction of the to- 
tal wave propagation time spent in the innermost regions, 
where growing modes are excited by the nonlinearity, scales 
with k2^n^ oc n~^, so the growth rate should scale also with 
rip^. Combining this with our observation that the growth 
rate increases approximately linearly with wave amplitude, 
we can write 

Im[cj](x^, (58) 

Tip 

for the modes excited when A < 1. 

Our simulations in BO 10 and Bll did not show any 
instabilities acting on the waves when A < 1. This can be 
neatly explained from Eq. 58 in the limit as ^ oo, where 
the growth rate tends to zero. This limit is appropriate since 
the waves have effectively no outer boundary in the simu- 
lations, because we damp the waves before they reach the 
boundaries of the computational domain. The fact that we 
observed no instabilities in the simulations is therefore con- 
sistent with this stability analysis, and is not a consequence 
of limited run time or insufficient spatial resolution. 



8.3 Comparison with the plane IGW problem 

This problem has some important differences with the case 
of a plane IGW in a uniform stratification. For that problem, 
as we discussed in the introduction to this paper, parametric 
instabilities act for any A, and always result in instability 
in the absence of diffusion. In addition, Lombard & Riley 
(1996) find that the presence or absence of isentropic over- 
turning does not seem to play a dominant role in, and is 
not the cause of, the instability^. This is different from our 
problem, where we find that overturning results in the pres- 
ence of a different class of localised modes, that are excited 
in the convectively unstable regions of the wave. However, 
we do find that the source of free energy driving the insta- 
bility, which is the free energy associated with primary wave 
entropy gradients, is the same whether A < 1 or A > 1. It 
is true that parametric instabilities exist for any A in our 
problem, like in the plane IGW problem, but these become 
unimportant in a large domain because the nonlinearity is 
spatially localised in the innermost wavelengths. This is dif- 
ferent from the plane IGW problem, in which the nonlinear- 
ity is important everywhere in the wave. 

The importance of overturning in our case could be be- 
cause the primary wave shear does not drive the modes, and 
in fact typically acts to stabilise them. In the plane IGW 
problem, instabilities for any A are driven by a combination 
of J^sw and Afbw , whereas in our problem instabilities are al- 
ways driven solely by J^tw Koudella Sz Staquet (2006) have 
also found jVbw to be the dominant energy source driving 
the instability of a (convectively stable) plane propagating 
IGW in 2D. They adopted a resonant triad model, which 
is valid for small wave amplitudes, to predict such a result. 



^ However, it is possible that in their calculations they have insuf- 
ficient resolution to be able to resolve any localised convectively 
unstable modes. If they go to larger A than the maximum they 
consider of 1.1, and/or consider larger resolutions, such localised 
modes may start to appear. 
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However, these calculations and ours were in two dimen- 
sions, so it remains to be seen whether shear would remain 
unimportant for our problem in 3D. 

From the results of their stability analyses, Lombard 
& Riley (1996) and Sonmor & Klaassen (1997) state that 
wave stability is a three-dimensional problem. This might 
suggest that the picture we have outlined could differ in 3D. 
However, the simulations performed in Bll show a strong 
similarity with the 2D results in BO 10. Performing a simi- 
lar stability analysis in 3D would be somewhat involved, and 
would be restricted to studying the stability of small ampli- 
tude waves, because the wave solution is not exact in 3D. 
It would be possible to calculate higher order terms to the 
solution, which would make it valid for larger A, and then 
perform a stability analysis of this wave. Without perform- 
ing such an analysis, it is difficult to quantify the importance 
of three-dimensional effects on the wave stability. Neverthe- 
less, the excellent correspondence between the results of the 
simulations in 2D and 3D suggest that for our problem the 
inclusion of a third dimension would be unimportant, with 
regards to wave stability. 



8.4 Implications for tidal dissipation 

In BIO and Bll, we discussed the implications of the wave 
breaking process for tidal dissipation in solar-type stars, and 
therefore to the survival of short-period planets. What more 
can we say in light of the results of this paper? 

One result is that parametric instabilities exist for waves 
with A < 1. These do not occur in an unbounded domain 
(the limit as rip oo), but will be present in the RZ of 
a solar-type star, since this does have an outer boundary, 
albeit many wavelengths from the centre of the star. We 
can roughly fit Eq. 58 to the results of our stability analysis, 
allowing us to give an upper bound to the expected growth 
rate of the strongest instability of a tidally excited gravity 
wave with A < 1. We write 



Im 



K- 



(59) 



and calculate a value of K from the solutions to our eigen- 
value problem, where we typically find K ^ 0.1. In the 
RZ of a solar-type star, tidally excited gravity waves have 
10^ < rip < 10^, for orbital periods in the range 1 ^ ^ 3 
days. We can therefore calculate an upper bound on the ex- 
pected growth rate of a parametric instability in a real star 
from taking A = 1 and Up - 10^, giving Im [uj/ujp] ^ 10~^, 
so that the resulting growth time. 



Im [cj] 



: 2.7 yrs. 



(60) 



It is important to note that this estimate is likely to be an 
approximate lower bound on tgrow , and will not be strongly 
affected by the inclusion of the rest of the RZ, because the 
the amplitude of the waves, and therefore the nonlinear- 
ity, is much smaller away from the centre. (In addition, our 
Boussinesq-type model is only valid where N (x r, which 
is only true near the centre of the star.) These calculations 
constrain the effects of nonlinear wave-wave interactions in 
the innermost regions, but do not take into account the rest 
of the RZ. 

It is important to estimate the magnitude of the re- 



sulting tidal dissipation, so that we can evaluate its role in 
the evolution of short-period planets. Instead of consider- 
ing the problem of continual forcing of the primary wave 
by the planet, we consider initialising the primary wave and 
ask how long it takes to be attenuated (and its energy dissi- 
pated), calling this timescale tni (this is similar to the highly 
eccentric binary problem discussed in KG96, which we dis- 
cuss in the next section). The torque on the star due to the 
gradual attenuation of the wave due to the combined action 
of these parametric instabilities at nonlinearly damping the 



(61) 



with the attenuation factor a = tgroup/tni, and F being 
computed as outlined in BIO and Bll. We define the global 
group travel time tgroup = 2 f^^ (l/cg^r)dr ~ 25 d, from a 
numerical calculation for the waves excited by a planet in 
a one-day orbit around the current Sun. The next question 
is: what is tnZ? To calculate this accurately is a very diffi- 
cult problem, and involves many uncertainties, particularly 
those involving the saturation process for these nonlinear 
couplings. However, we note that a lower bound on tni can 
be obtained by the growth rate of the fastest growing para- 
metric instability tgrow This is because this will act as a 
bottleneck for the nonlinear cascade of energy from the pri- 
mary wave, and so will limit the maximum decay rate of the 
primary wave. This is probably also true if we are continu- 
ally forcing the wave. We can then estimate 

a < ^ 0.025. (62) 

tgrow 

This gives an upper bound on the torque resulting from the 
nonlinear damping of the primary wave. 

The resulting tidal quality factor can be computed from 
the expression 



4T V -h rrir 



m^Rl /^27ry 



dyn 



P ' 



(63) 



where m^,p are the stellar and planetary masses, is the 
stellar radius and oo'^yn is the square of the dynamical fre- 
quency of the star Grrii./R^. As before, P is the planetary 
orbital period. This can be used to give a lower bound on 
the tidal quality factor resulting from nonlinear damping of 
the primary wave in the A <1 regime, where we find 

^5 



10^ 



1 



' a"^10^ ^bx 10^ 



(64) 



in the weak damping limit. The efficiency of this process 
is less than critical layer absorption by a factor ^ 1. 
Note that this gives a lower bound on Q'^^ because tni is 
likely to be somewhat larger than tgrow (e.g. KG96 take 
tni = ^Otgrow)- In addition, we have taken the most opti- 
mistic value of A = 1, corresponding to the waves excited by 
a planet with a mass of about 3Mj (see Eq. 1). The result- 
ing may therefore be one or several orders of magnitude 
larger than this lower bound. Furthermore, it is interesting 
to note that this bound may not be sensitive to the num- 
ber of wavelengths in the RZ, and therefore to the orbital 
period, because the dependences of tgroup and tgrow on rip 
cancel at leading order. 

The parametric instabilities that exist when A < 1 are 
much slower than the rapid instabilities that onset when 
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A > 1. The nonlinear outcome of the A > 1 instabihties is 
that the wave breaks and forms a critical layer, which then 
absorbs subsequent ingoing waves, and results in astrophys- 
ically efficient tidal dissipation. The estimate of this section 
indicates that the parametric instabilities that exist when 
A < 1 are much less efficient at dissipating energy in the 
tide, by several orders of magnitude. This important result 
supports the explanation outlined in BIO and Bll for the 
survival of short-period planets around solar-type stars. 

8.5 Comparison with Kumar Goodman 

We can qualitatively compare our results with previous work 
by KG96, who studied nonlinear damping of tidal oscilla- 
tions in highly eccentric solar-type binaries. They used a 
truncated Hamiltonian approach to study parametric insta- 
bilities of tidally excited f and g-modes. In their model, stel- 
lar eigenmodes are coupled together from terms that ex- 
ist at third order in displacement in the expansion of the 
Lagrangian density, i.e., they adopt a weakly nonlinear ap- 
proach. They consider the evolution of a mode that has been 
tidally excited, but is no longer subject to forcing, due to 
nonlinear coupling with a large number of g-modes that are 
present in the RZ of the star (and exist because they have 
already been excited by turbulent convection, for example). 
Their result indicates that high order and high degree g- 
modes can be parametrically excited by low order quadrupo- 
lar f and g-modes, and can draw energy from the primary 
mode on a timescale that is much shorter than the radiative 
damping time of the primary mode. 

A direct comparison of our work with theirs is not possi- 
ble for several reasons. Firstly, in their numerical work they 
mainly consider a primary f-mode coupling to many g-modes 
in the RZ. The f-mode eigenf unction has its largest magni- 
tude at the surface and decays rapidly inwards, in contrast 
to the primary g-modes that we are considering, so the cou- 
pling strengths are likely to be different. Secondly, we only 
consider the nonlinear interactions in the central regions of 
the star, where they are likely to be most important for 
g-modes, whereas they consider these interactions through- 
out the whole star. Thirdly, our model is 2D, whereas their 
eigenfunctions are valid in 3D for a spherically symmetric 
background. This last point, however, is probably not im- 
portant. 

One important point is that they neglect the possibility 
of wave breaking, which would provide an upper limit to the 
amplitude of a given mode. This would prevent modes with 
large amplitudes from coupling with the primary wave, and 
the nonlinear outcome of the breaking (most notably critical 
layer formation) would significantly modify the strength of 
tidal dissipation. Their results will therefore not be valid for 
primary or daughter waves that satisfy a breaking criterion, 
since weakly nonlinear theory is insufficient in this case. In- 
deed, the concept of parametric instability is no longer valid 
if the daughters break and cannot form standing modes. 
This is particularly important given their primary applica- 
tion of eccentric solar-type binaries, since in that case, the 
ampliudes of the waves are likely to be large enough for wave 
breaking near the centre of their stars (this is estimated in 
the Appendix of OL07, for example). 

Keeping in mind the differences between our approach 
and theirs, we now directly apply their results to our prob- 



lem, and quantitatively compare the growth time of para- 
metric instabilities with those found in this paper. The 
growth time in their work 



V 1035 J J 



-1/2 



(65) 



where Ep^o is the initial energy in the primary wave. For the 
g-modes that we consider. 



'p,0 



J J J Sr'^ sin OdrdOdcj) 



(66) 



This can be computed to give Ep^o ^ 2x 10^^ J for a Jupiter- 
mass planet on a one-day orbit around the current Sun, 
which has A ^ 0.3. This means that tgrow ~ 3 yr when 
A = 1, which happens to agree surprisingly well with our 
calculation in the previous section, given the differences in 
our approach. The total number of daughter modes which 
simultaneously interact with the primary in their model is 

- 10^° (^^^y^^ ^ 10^ for our fiducial case. We also find 
that there are many growing modes for a given set of param- 
eters in our stability analysis, so these statements appear 
qualitatively consistent. They find that collectively, these 
modes absorb most of the energy of the primary wave after a 
time ~ lOtgrow (this is equivalent to assuming a ~ 2.5xl0~^ 
in the previous section). This predicts Q^^ 5 x 10^ in their 
approach. We therefore conclude that our results are broadly 
consistent with KG96. 



9 CONCLUSIONS 

In this paper we have performed a stability analysis of a 
standing internal gravity wave near the centre of a solar- type 
star, using the 2D exact wave solution derived in BO 10. This 
work has relevance to the tidal interaction between short- 
period planets and their solar-type host stars, since these 
waves are excited at the top of the radiation zone of such a 
star by the tidal forcing of the planet. The equations gov- 
erning the evolution of the perturbations to this wave were 
written down in spectral space using a Galerkin spectral 
method, and then solved as an eigenvalue problem. This 
required the imposition of an artificial impermeable outer 
boundary several wavelengths from the centre of the star. 

We have identified the modes that initiate the break- 
ing process when the wave overturns the stratification. This 
type of mode is strongly localised in the convectively unsta- 
ble regions of the primary wave, and is driven by unstable 
entropy gradients. Its growth time is comparable with the 
primary wave period, which is consistent with the breaking 
time observed in the simulations of BO 10 and Bll. 

We have also studied the instabilities which exist for 
waves with insufficient amplitudes to overturn the stratifica- 
tion. We find that these are parametric instabilities driven 
by (convectively stable) entropy gradients in the primary 
wave. The growth rate of these modes scales inversely with 
the number of wavelengths within the domain, so they be- 
come less important for a real star than for the small con- 
tainer considered here. It is estimated that their growth 
times in a real star would be of the order of 3 yr, which is 
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much longer than the orbital period of a short-period planet, 
though many such modes are excited. Rough estimates are 
made that provide an upper bound on the resulting tidal 
dissipation, for which we find a lower bound on the tidal 
quality factor of Q'^ > 10^ from this process. This is much 
weaker than the dissipation resulting from critical layer ab- 
sorption obtained in BO 10 and Bll, and so is unlikely to 
change the picture outlined in Bll for the survival of short- 
period planets. 

The results of this paper provide further support for 
the hypothesis outlined in BO 10 and Bll for the survival of 
short-period extrasolar planets around slowly rotating solar- 
type main-sequence stars. Coupled with weak dissipation of 
the stellar equilibrium tide by turbulent convection when 
the orbital period is shorter than the convective timescale 
(see e.g. Zahn 1966; Goldreich & Nicholson 1977; Goodman 
& Oh 1997; Penev & Sasselov 2011), and the absence of 
inertial wave excitation in their slowly rotating stellar hosts 
(OL07), it seems likely that short-period planets can survive 
against tidally induced orbital decay if they are unable to 
cause the internal gravity waves that they excite to break 
near the centre. This paper demonstrates that the waves 
need to overturn the stratification near the centre to obtain 
efficient tidal dissipation. 

We discussed several differences between our problem 
and the stability of a plane IGW in a uniform stratiffcation 
(e.g. Lombard & Riley 1996). We have confirmed that when 
the wave is confined in a container with an outer boundary 
it is unstable whatever its amplitude, in the absence of dif- 
fusion. However, the inverse dependence of the growth rate 
on the number of wavelengths within the container is quite 
different, and results from the finite time of nonlinear inter- 
action being much shorter than the group travel time across 
a large container. 

We compared our results to Kumar & Goodman (1996), 
who studied the nonlinear damping of tidally excited oscilla- 
tions in highly eccentric binaries, and found some agreement. 
They predict that many (~ 10^) modes collectively draw 
energy from the primary wave, which we have qualitatively 
confirmed from our stability analysis. The growth rates of 
parametric instabilities for the same problem in both of our 
approaches when A ^ 1 are very similar. They therefore pre- 
dict a similar lower bound for resulting from this process. 
This is promising, given the differences in our approach. It 
would be interesting to extend their numerical calculations 
by studying the parametric instabilities of g-modes includ- 
ing continual tidal forcing of the primary wave and nonlin- 
ear couplings involving many daughter and granddaughter 
modes, as well as taking into account the amplitude limit- 
ing effects of wave breaking. Weakly nonlinear theories such 
as ours and theirs are likely to be valid when considering 
the initial stages of the breaking process, and in study- 
ing whether any instabilities exist for suboverturning waves, 
which were the topics of study in this paper. However, they 
should not be used to determine long-term behaviour for 
waves which overturn the stratification (whenever A > 1). 
This means that for the circularisation of eccentric solar- 
type close binary stars, it is inappropriate to use a weakly 
nonlinear approach, since in that case wave breaking is very 
likely to occur. Instead, the results of BO 10 and Bll must 
be used to obtain the correct magnitude of the dissipation, 



and the resulting circularisation rate due to nonlinear inter- 
actions between gravity waves. 

It would be worthwhile to confirm the results of this 
paper using 2D numerical simulations with SNOOPY, such 
as those described in BO 10. An artificial impermeable outer 
boundary could be implemented in the code, and the result- 
ing instabilities then studied. Of particular importance is to 
determine the rate at which energy is lost from the primary 
wave due to the parametric instabilities for suboverturning 
waves that we studied in this paper (i.e., to numerically cal- 
culate tni)- This would enable a more accurate calculation 
of the magnitude of Q'^ and would provide a useful inde- 
pendent check of our results. We defer such calculations to 
future work. 
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APPENDIX A: TOY MODEL: PARAMETRIC 
INSTABILITY OF PRIMARY WAVE 

Parametric instability is a type of resonant triad interaction 
in which the transfer of energy from a parent (subscript p) 
mode, with amplitude Ap, destabilises a pair of daughter 
(subscript dl,d2) modes (which exist when Ap = 0). These 
can then be damped or subject to further nonlinear inter- 
actions (to produce granddaughter modes, and so on). The 
frequencies of the modes must satisfy an approximate tem- 
poral resonance condition ujp ^ ujdi + ujd2 , for parametric 
resonance to occur. 

The equations governing the temporal evolution of the 
mode amplitudes take the form (e.g. Dziembowski 1982; Wu 



& Goldreich 2001) 

Ap = ^pAp — iujpAp -\- iujpaAdiAd2, (Al) 

Adi = -jdiAdi - iudiAdi + iudicrApAd2, (A2) 

Ad2 = -'yd2Ad2 - ioJd2Ad2 + iuJd2crAdiAp. (A3) 



In these equations, 7^ is the linear growth/damping rate of 
mode J, and cr is the nonlinear coupling strength for these 
three modes. Here Aj is the amplitude of mode j, with the 
energy in that mode being proportional to \Aj\^ . 

The coupling coefficient a is largest when Udi ~ 00d2 = 
uj, and therefore uj ^ ujp/2. li the daughter modes have sim- 
ilar frequency, then we can assume that they have similar 
spatial scales. Hence we can take their damping rates to be 
the same, i.e., 7^1 = 7^2 = 7- To consider the initial stages 
of the breaking process we take Ap to be approximately con- 
stant in time. In our problem the primary wave is maintained 
at a constant amplitude due to forcing, and is not unstable. 
The evolutionary equations reduce to 

Adi = -^Adi - iojAdi + i(jjaApAd2, (A4) 
Ad2 = -'^Ad2 - iojAd2 + iojaAdiAp. (A5) 
If we take Adi oc expst, then the growth rate is 

Re[s] = -^ + ^cva\A,\. (A6) 
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The 



growth rate is reduced if the detuning Au = 
^di — ^d2 / 0, by changing the second term to 



{l/2)^^UJdlUJd2cr'^\Arp\'^ — (Acl;)^. From this model, we expect 
7 to simply reduce the growth rate for a given mode. In 
addition the growth rate scales linearly with the amplitude 
of the primary (parent) mode. The threshold amplitude for 
instability in this simple model is 



7 



(A7) 



which depends on the coupling strength cr. 

The spatial dependence of the interaction is contained 
in the coupling coefficient a, which contains an integral of 
the product of the three eigenf unctions. This toy model of 
parametric instability is useful as a simple model to under- 
stand some of the results of § 6. It is interesting to note that 
in this model, uj ^ m/km^n ^ for A ^ 1, since in this 
limit the daughter modes have frequencies comparable with 
the linear mode frequencies. This results in a growth rate 
scaling inversely with Tip. 



APPENDIX B: A SINGLE IGW IS IN 
EQUIPARTITION 

An IGW with a single value of m and n satisfies equiparti- 
tion of kinetic and potential energies, when integrated over a 
multiple of half- wavelengths, as we will now prove. If we take 
f{r) = Jm{km,nT)e^^^ ^ wc cau rewrite Bessel's equation in 
the form 

(Bl) 



1 



After multiplying by r/, and then integrating over radius 
from ri to r2, we obtain 



•^ri L 



rdr 



= r kl,^nfdr+[rfdrfYrl.{B2) 
J n 



Since 
K 



J n 



2 



rdr, 



is the integrated kinetic energy, and 



rr2 ^ 

J n 



(B3) 



(B4) 



is the integrated potential energy, as defined in the text, 
this statement is telling us that equipartition holds if we 
integrate over a range where / or drf are zero at the end 
points. 
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